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Abstract Partially implicit Runge-Kutta methods are presented in this work in 
order to numerically evolve in time a set of partial differential equations. These 
methods are designed to overcome numerical instabilities appearing during the 
evolution of a system of equations due to potential numerical unstable terms in 
^ 1 the sources, such as stiff terms or the presence of factors as a result of a partic- 

<n: ular chosen system of coordinates. In this article, partially implicit Runge-Kutta 

methods for several convergence orders have been derived and stability properties 
have been analyzed. These methods are shown to be appropriated to avoid the de- 
CL, velopment of numerical instabilities in the evolution in time of wave-like equations 

Jh \ in spherical-type coordinates, in contrast to the explicit Runge-Kutta methods. 
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<^^> 1 Introduction 

CO 
Ov 

The evolution in time of many complex systems, governed by partial differential 
equations (PDE), implies, in a broad variety of cases, looking for the numerical 
solution of a system of ordinary differential equations (ODEs). The most commonly 
used methods to integrate in time these systems of ODEs are the well-known 
Runge-Kutta (RK) ones (see e.g. [1] for a general review of these methods and their 
main properties) . Several classifications of the RK methods can be done, according 
to, e.g., their convergence order, the number of stages or their explicit /implicit 
structure. 
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Implicit or partially implicit methods are used to deal with systems which 
require a special numerical treatment in order to have a stable evolution. The 
origin of the numerical instabilities in the evolution of a system of equations is 
diverse. Stiff terms in the sources of the equations can lead to the development of 
numerical instabilities when explicit RK (ERK) methods are used. In this context, 
the so-called implicit-explicit (IMEX) RK methods have been used to evolve in 
time conservation laws with stiff terms or convection-diffusion-reaction equations 



Not only stiff terms can be the reason of the appearance of numerical insta- 
bilities. Source terms of the equations, when a particular system of coordinates is 
chosen, may introduce factors which can also be numerically interpreted as stiff 
terms. This fact will be evident for a wave equation in spherical coordinates in 
the numerical examples shown in Sect. [H in which the l/r p factors, p £ N, ap- 
pearing due to spherical coordinates are interpreted as stiff terms close to r = 0, 
even when the evolved data is regular. In the case of the development of these 
numerical instabilities, an implicit treatment of the system offers a solution to get 
a stable evolution. 

The implicit treatment of the sources needs, in general, an inversion of some 
operators. Depending on the complexity of the equations, the inversion can be done 
analytically or numerically, or be even prohibitive in practice from the numerical 
point of view. We will focus on a particular structure of equations which does 
not require any analytical or numerical inversion. Therefore, these methods have 
a computational cost similar to the ERK methods but they are able to provide 
stable evolutions due to their partially implicit component. 

The choice of the coefficients in the derivation of the methods presented in this 
work are based on stability properties for both explicit and implicit parts, as it is 
described in the following sections. The stable evolution obtained for a particular 
example shows their validity. However, potential applications in other systems of 
equations have still to be checked. 

The manuscript is organized as follows. In Sect.[2]the structure of the system of 
equations and the requirements for the numerical scheme are described. In Sect.|3l 
partially implicit RK (PIRK) methods, from first to third order, are derived; some 
comments concerning the derived PIRK methods and some available IMEX ones 
are mentioned. Numerical experiments are shown in Sect.[H where PIRK methods 
are used to evolve in time a wave equation in spherical coordinates. Conclusions 
are drawn in Sect. [5j 

2 Structure of the equations and requirements for the numerical methods 

Let us consider the following system of PDEs, 



being C\, £2 and £3 general non-linear differential operators. Let us denote by L±, 
L2 and L3 their discrete operators, respectively. No restrictions onto the discrete 
operators are imposed. They can correspond, for example, to a flux-conservative 
numerical scheme of a conservation law or a finite difference scheme of a general 
evolution system. 




(i) 
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L\ and L3 will be treated into an explicit way, whereas the L2 operator will be 
considered to contain the unstable terms and, therefore, treated partially implic- 
itly. Each stage of the RK method will proceed into two steps: i) the variable u 
is evolved explicitly; ii) the variable v is evolved taking into account the updated 
value of u for the evaluation of the L2 operator. This strategy implies that the 
computational costs of the methods are comparable to those of the explicit ones. 
The resulting numerical schemes do not need any analytical or numerical inver- 
sion, in the same way as ERK methods do not, but the partially implicit treatment 
confers stability to the numerical algorithm. 

Numerical methods based on a nonlinear stability requirement are very desir- 
able. Such methods were originally named as total variation diminishing (TVD) 
and are also referred to as strong stability preserving (SSP) methods (see e.g. [6j 
[7]). If U = U(t) is a vector of discretized variables, i.e., [U(t)]j = Uj(t) = u(xj,t), 
and it™ is the numerical approximation to u(xj,t n ), then SSP discretizations have 
the property that the total variation 

TV(U n ) = J2W-^-i\ (2) 
3 

of the numerical solution does not increase with time, i.e., 

TV(U n+1 ) < TV{U n ). (3) 

A sequence U n is said to be strongly stable in a given norm ||<|| provided that 
\\U n+1 \\ < || (7™ || for all n > 0. The choice of the norm is arbitrary, being the 
TV-norm, and the infinity norm, two natural possibilities. 

A s-stage ERK method for the equation dtU = L(U) can be written in the 
form: 

t/ (i) = l ^(a ik U^ + Atp ik L{U^)) , 

fc=0 

i = 1,2, . . . ,s, 

U n+1 =U {s \ (4) 

where all the > 0, and = only if /3j^ = 0, and At denotes the time step. 
Gottlieb and Shu !6] proved that the classical second-order method, 

f/(°) = W n , 

v (i) = u n + AtL(U n ), 

U n+1 = i[/" + ic/ (1) + iz\tL(t/ (1) ), (5) 

is the optimal second-order two-stage SSP ERK method, and that the third-order 
one due to Shu and Osher [8], 

f/ (0) = U n , 

f/ (1) = U n + AtL(U n ), 

UW = lu» + luW + lAtL{uM), 
4 4 4 v 

U n+1 = ~U n + \u { - 2) + ^AtL(U (2) ), (6) 
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is the optimal third-order three-stage SSP ERK method. The optimal adjective 
refers, for a given number of stages, to a maximization of the corresponding 
Courant-Friedrichs-Lewy (CFL) value (1 in both cases) and the efficiency in the 
storage requirement. 

This property is crucial to capture the right behavior of the evolution of the 
system. In the derivation of the PIRK methods proposed in this work, the previ- 
ously described optimal SSP RK methods are recovered when the L2 operator is 
neglected, i.e., when partially implicitly treated parts are not taken into account. 

Positiveness of the coefficients is preferable. See [IJ[9] for more details about the 
implications of positive/negative coefficients on stability properties and storage 
requirements. The PIRK methods derived minimize the number of stages, two 
and three in the case of second and third-order methods, respectively. They differ, 
among other aspects, from the explicit, singly diagonally implicit RK ones [10] by 
not having an explicit first stage. Taking into account previous considerations, the 
remaining coefficients associated to the L2 operator for the different methods are 
chosen according to stability criteria. 

The requirements for the numerical algorithm to solve the system ([TJ) are sum- 
marized in the following points: 

i) RK method, with at most a singly diagonally implicit RK scheme for the 
implicit part. 

ii) Positive coefficients are preferable. 

iii) Recovery of the optimal SSP ERK methods when implicitly treated parts are 
neglected. 

iv) Choice of the remaining coefficients according to stability criteria. 
3 Numerical methods and stability analysis 

In this section, the PIRK methods, from first to third order, are derived. The meth- 
ods for a particular convergence order and equal number of stages, and imposing 
SSP optimal ERK methods for the pure explicitly treated parts are deduced; next, 
stability analysis is carried out and optimal values for the remaining coefficients 
are selected, with preference to positive coefficients. 

3.1 Stability definitions 

Numerical methods combining explicit and implicit treatments can be analyzed in 
several ways from the stability point of view. Let us consider the scalar complex 
test equation, 

d t y{t) = uy{t) + f iy(t), (7) 

where y(t) is the evolved variable, and v and fl denote the eigenvalues correspond- 
ing to the explicit and implicit operators, respectively. Let us define v :=v At and 
fi := flAt, where At is the time step of the RK method. When a particular RK 
scheme is applied to Eq. (J7J), the updated value of the evolved variable, y n , can 
be written in terms of its previous value, y n , as 



n+l r,/ \ 

y = R{n,v)y 



(8) 
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where R(n, v) is a fraction of polynomials depending on fj, and v, called stability 
function. Stability thus requires that \R(fi,u)\ < 1. 
Let us define 

T>q : = G C : the scheme is stable for any v g C~ j . (9) 

Here, stability with respect to the implicit part is considered, and it has to be 
compared whether T>q is smaller than the stability region of the explicit method. 
Alternatively, one can insist on using the full stability region of the explicit method, 
and, in that case, the following set has to be taken into account: 

Vi := | J/ G C : the scheme is stable for any fi G C-,\l+n\<l\. (10) 

In this article we will focus on the set T>\ from now on. A SSP ERK method 
is recovered when implicitly treated parts are neglected. This property guarantees 
the stability of the pure explicit parts. The stability analysis will concern the 
addition of the partially implicit component. 

Since we are considering a system of equations, and the characteristic structure 
of eigenvalues and eigenvectors of the explicit and implicit parts do not coincide 
necessarily, the previous analysis has to involve matrices instead of scalars, and 
the global structure of the system has to be taken into account. Let us denote 
(cviu, a2v), Am and (71M, 72^) the associated linearized parts of the C\, £2 and £3 
operators, respectively. The linearized system ([I]) is then rewritten as 

u t = a%u + a 2 v, , , 

vt = 7111 + 72^ + Ait. 

Let us denote a, := ctiAt, A := XAt and 7, := jiAt. A stability function and its 
absolute value will be replaced by a stability matrix and the absolute values of 
its eigenvalues. The hypothesis [1 + fj,\ < 1 will be replaced by [tOj| < 1, where w,, 
i = l,2 denote the two eigenvalues of the matrix regarding the explicit part 



1 + aj 02 
71 1 + 72 



(12) 



This condition implies, in particular, that the determinant of previous matrix, 
denoted by dex, and its trace, denoted by trex, are bounded as follows: 

|dex| = |wi||w2| = |(1 + qi)(1 +72) - «27i| < 1) (13) 
|trex| = \uj\ + lj 2 \ = |2 + ai + 72 1 < 2. (14) 

Moreover, in the condition \R(n, v)\ < 1, the stability function R(fi,u) will be 
replaced by the two eigenvalues associated to the matrix M, M defined such that 

""' ' ' l m( «" ), (I.:,) 



V 



for a given method. However, in order to simplify the derivation of the PIRK 
methods, we are going to relax the bound on the eigenvalues of M by a bound 
on its determinant, i.e., | det(M)| < 1. The restriction onto the eigenvalues will be 
recovered in the numerical experiments in Sect. [3] as the boundaries of the stability 
region. 

We also assume that Re(Aa2) < 0, where Re denotes the real part of a com- 
plex number. This condition is satisfied for general sinusoidal wave-like equations 
written as a first-order system in time (see Sect. H] for examples). 
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3.2 First-order method 

The derived one-stage first-order method for the system (Q} can be written in 
terms of one coefficient, ci, as follows: 



U =U + At L\(U ,V ), 

n+1 
V = V 



n + At [(1 - Cl )L 2 {u n ) + Cl L 2 {u n+1 ) + L 3 (u n , v n )] 



(16) 



A Taylor expansion of (u n+1 , « n+1 ) and all the operators in terms of (u n ,v n ) and 
their spatial derivatives was carried out to guarantee the convergence order of 
the method. This method is a particular case for the system ([1]) of the so-called 
IMEX-0 method (see, e.g., [TT]), where c± is the 6 parameter. 
System (|16[) can be written as 



v 



Mr I " I, (17) 



where 



Mi=[ In" 1 v , a \ I (18) 

' 71 + A(l + aici) 1+72 + Aq 2 ci 



Its determinant is: 



det(Mi) = dex- Aq 2 (1 - ci). (19) 

Taking into account bound (|13|) . c± = 1 guarantees | det(Afi)| < 1, for all possible 
values of Aa 2 . The resulting method is: 



u =u + At L\{u ,v ), 



71+1 

V 



v n + At [L 2 (u n+1 ) + L 3 (u n ,v n )] 



(20) 



3.3 Second-order method 



The derived two-stages second-order method for the system ([1), imposing SSP 
optimal two-stages second-order method for the pure explicit parts (LI and L3 
operators), can be written in terms of two coefficients, {c\,C2) (do not confuse 
coefficient c\ with the one appearing in the first-order method), as follows: 



u w =u n + AtL 1 {u n ,v n ) 



71+1 



n+1 



v n + At [(1 - Cl )L 2 {u n ) + ciL 2 (« (1) ) + L 3 (u n ,v n )] 

:^[u n + u^+AtL 1 (u^,v^)], 
v n + ^ [L 2 (u n ) + 2c 2 L 2 (u^) + (1 - 2c 2 )L 2 (i 
+L 3 (u n ,v n ) + L 3 (u^,v ( - 1) )]. 



"+1\ 



(21) 



(22) 



A Taylor expansion of (u^\ v^ 1 '), (u n+1 , v n+1 ) and all the operators in terms of 
(u n , v n ) and their spatial derivatives was carried out to guarantee the convergence 
order of the method. 
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From (ED and (EZ 



,n+l 
,,n+l 



Mo 



(23) 



where 
M 2 = 



1 

(l/2-c 2 )A 1 



1 



+ - 



1 + Q?l Ct 2 

71 + 2Ac 2 72 



1/2 
(A + 7 l)/2 1 + 72/2 

1 + a\ a 2 
7i + A(l + cnci) 1 + 72 + Aa 2 ci 



(24) 



Its determinant is: 
det(M 2 ) = 



i [(1 - dex) 2 + trex 2 + Aa 2 (l - dex)(l - 2ci + 2c 2 ) 
+A 2 a|(2c 2 - ci - 2cic 2 ) . 



(25) 



The conditions 1 — 2ci + 2c2 = 2c2 — ci — 2ciC2 = guarantee | det(M2)| < 1, for all 
possible values of \a 2 . But they lead to ci = (1 ±i)/2 ^ R, C2 = ±j/2 ^ R, where 
i = is the imaginary unit. 

Let us restrict to real numbers, so coi £ R and Act2 £ R~. j det(M2)j < 1 is then 
equivalent to: 



where 
and 



4 < Ki + Xa 2 K 2 (l - lex + 2c 2 ) + A 2 al(2c 2 - ci - 2cic 2 ) < 4, 



Ki := (1 - dex) 2 + trex 2 = 1 + u)\u)l + u\ + w| G [1,4] 



(26) 



(27) 



if 2 := 1-dexe [0,2]. 



(28) 

The minimum value of K\ is reached for ui\ = uj 2 = 0. The maximum value of K\ is 
reached for ui\ = ui\ = 1. The minimum value of K 2 is reached for uj\ = ui 2 = ±1. 
The maximum value of K 2 is reached for uj\ = —uj 2 = ±1. For ui\ = uj 2 = 0, 
K 2 = 1. 

We analyze the value of det(M2) for cjj = 0, ±1, i = 1,2, and in the cases 
I Acv2 1 S> 1, Aa?2 ~ — 1 and \Xa 2 \ <C 1. The resulting sufficient conditions are (see 
Apendix [Al for more details): 



together with 



< 2c 2 (l -a) -a, 1 - 2ci + 2c 2 < 0, 
< 6 + 5ci - 6c 2 + 2cic 2 , < 4 + ci - 2cic 2 , 



4 < Aq 2 (1 - 2ci + 2c 2 ), -5 < A 2 a^(2c 2 - ci - 2cic 2 ). 



(29) 



(30) 



For |A«2| <§C 1, first inequality in (|30p is the most relevant one of both. We 
impose 1 — 2ci + 2c2 = , and minimize |2c 2 (1 — ci ) — ci | (taking into account t|29|) ) . 
The resulting values are c\ = 1/2 and c 2 = 0. The scheme is then written as: 



u (1) =u n + AtL 1 (u n ,v n ) 
(1) = v n + At 



V 



i L2 (u n ) + ±L 2 (u^)+L 3 (u n ,v r - 



(31) 
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n+1 



n+1 



1 \ u n + u^ +AtL 1 (u^,v^ 

A 

77. ^1 t / 7J. \ t" / TJ.-M 



'* i 



L 2 (u n ) + L 2 (u n+1 ) + L 3 (u n ,v n ) + L 3 (« (1) ,« (1) )] 



(32) 



For |A«2| 3> 1, second inequality in (|30[) is the most relevant one of both. We 
impose 2c2 — ci — 2ci C2 = 0, and minimize 1 1 — 2ci + 2c2 1 (taking into account (|29|) ) . 
The resulting values are c\ = 1 — v2/2 and C2 = (\/2 — l)/2. The scheme is written 
in this case as: 



n+1 



" + 4t^(«",/), 
V2 



v n + At 



L 2 (u n )+[ 1 



V2 



2 



J L2^ (1) ) + i 3 K,«") 



i L" + u (1) +ZkLi(u (1) ,« 



=«" + ^ |l 2 (u") + (^-1)X 2 (u«) + (2-V2> 2 (u" +1 ) 



(33) 



(34) 



+L 3 (n", U n )+L 3 ^ (1) ^ (1) )] 



Depending of the value of |Aa2| = Aa^l^ 2 , it will be more convenient to use 
a particular set of values for the coefficients. For not too stiff numerical problems, 
the choice (c\,c 2 ) = (1/2,0) will be more appropriated and it avoids to compute 
L 2 (vS 1 ') to obtain v n+1 in the final stage. For very stiff numerical problems, last 
option is better. This fact will be illustrated in Sect. [3] 



3.4 Third-order method 

The derived three-stages third-order method for the system (fTJ), imposing SSP 
optimal three-stages third-order method for the pure explicit parts (LI and L3 
operators), can be written in terms of two coefficients, (c\,c 2 ) (do not confuse 
these coefficients with the ones appearing in previous methods), as follows: 

u (1) =u n + AtL 1 (u n ,v n ), 

=v n + At [(l- Cl )L 2 (u n )+ Cl L 2 (u^)+L 3 (u n ,v n )} , (35) 

' = l^3u n + u^ +AtL 1 (u^,v^)], 

< „(3) = v n + A* r 2(ci + 2c2 )L 2 {u n ) + 4c2L 2 (m (1) ) + 2(1 - ci - 4 C2 )L 2 ^ (2) ) 
+L 3 (u n ,v n ) + L 3 (u ( - 1 \v^)], 

(36) 

' u n+i = ~[u n + 2u^ + 2AtL 1 (u^,v^)] , 

« v n+1 =v n + ^[L 2 (u n ) + L 2 (u ( - 1 '>) + 4L 2 (u^) (37) 
+L 3 (u n ,v n ) + L 3 (u ( - 1 \v^)+4L 3 (u < - 2 \v^ ) )], 

A Taylor expansion of the values of the variables u and v in all the stages of the 
method and all the operators in terms of (it™, v n ) and their spatial derivatives was 
carried out to guarantee the convergence order of the method. 



Partially implicit Rungc-Kutta methods 



9 



From (I35l)-(l3 



where 



u n+1 



M 3 ( ..„ , (38) 



N 2 



(l-ci-4c 2 )| 1 



1 + ^ 



4 



^ + ( Cl +2c 2 )^ 1 + f 



+ 



ILL 

4 



iii 

4 



(41) 



Its determinant is: 

det(M 3 ) = — [l4 + 2(trex - l) 3 + (dex - 2) 3 + 6 trex 2 + 3 dex((trex - l) 2 - 2) j 
+ 2^A« 2 (-1 + ci - 4c 2 ) [(dex - 2) 2 + (trex - l) 2 - 2 
+ ^A 2 a| fci - 4c 2 + (dex - l)(4c 2 - c 2 - 4cic 2 
1 rA 3 a 2 [-l + 3(l-2ci)(ci+4c 2 )]. 



72 



(42) 



The expressions — l+c±— 4c 2 , ci— 4c 2 , 4c 2 — cf — 4cic 2 and —1+3(1— 2ci)(ci+4c 2 ) 
cannot vanish simultaneously. Let us restrict, as in previous subsection, to real 
numbers, so uji G ffi and Aa 2 g R~. | det(Ms)| < 1 is equivalent to: 

_ K3 Aq 2 ^ 4 (-1 + ci -4c 2 ) 
- 36 24 



A OH 

+^2^[ci - 4c 2 + (dex - l)(4c 2 - c? - 4cic 2 ))] 

\ 3 _ 3 

-l + 3(l-2ci)(ci + 4ca)] < 1, 



72 



(43) 



where 

K 3 := 14 + 2(trex-l) 3 + (dex-2) 3 +6trex 2 + 3dex[(trex-l) 2 -2] G [-12,36], (44) 
and 



K 4 : = (dex - 2) 2 + (trex — l) 2 — 2 G [0, 



(45) 



The minimum value of K3 is reached for wi = 1 = — oj 2 and 102 = 1 = — wj.. For 
cji = oj 2 = — 1, K3 = 4. The maximum value of K3 is reached for wi = w 2 = 1. 
The minimum value of K4 is reached for uj± = cj 2 = 1. The maximum value of K4 
is reached for uj\ = w 2 = — 1, wi = 1 = — w 2 and oj 2 = 1 = — cj%. 
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We analyze previous bounds for the values uii = ±1, i = 1,2, and in the cases 
|Aa 2 | 3> 1, Aa 2 w — 1 and |Aa 2 | <S 1. The resulting sufficient conditions are (see 
Apendix [B] for more details): 

20 

_ — < Cl - 4c 2 < 0, -1 + 3(1 - 2ci)(ci + 4c 2 ) < 0, 

< 73 + 18c? - 180c 2 + 9ci(3 + 8c 2 ), < 9ci - 12c 2 - 6c? - 24cic 2 + 143, 
< 103 - 15ci - 6c? + 84c 2 - 24cic 2 , < 6c? - 15ci + 36c 2 + 24cic 2 + 71, 



(46) 



together with 



Aa 2 (-1 + ci - 4c 2 ) < 



-24 < A 2 a|(ci - 4c 2 ), 



A 3 o|[-l + 3(1 - 2ci)(ci + 4c 2 )] < 48. 



(47) 



For A« 2 | <C 1, first inequality in (|47|) is the most relevant one. Firstly, we 
minimize | — 1 + ci — 4c 2 | (taking into account (|46p ). Consequently, we choose c 2 = 
ci/4. This condition minimizes the factors accompanying Aq 2 and A 2 q 2 in (|47|) . 
The remaining inequalities from (|46[) and (|47(1 reduce to 



3 - ^ < Cl < 3 + ^, A 3 a|(-1 + 6 C1 - 12c?) < 48. (48) 



12 



12 



Secondly, we minimize | — 1 + 6ci — 12c?j (taking into account its allowed rang 
The minimum is placed at ci = 1/4. Therefore, the resulting values are 



(ci,c 2 ) 



1 _L_ 

4' 16 



(49) 



and the method is written as 



■ U -\- At L\{U ,V 

v n + At 



\L 2 {u n ) + \L 2 {u^)+L a {u n ,v n ) 



(50) 



(2) 



„ . At 

v +T 



|l 2 ( w ») + Il 2 ( w W) +j L 2 ( u ( 2 )) 
"+L 3 K,^)+i3(« (1) ,« (1) )] , 



(51) 



,n+l 



n+1 



1 \u n + 2u^ +2AtL 1 (u^,v^) 
«" + ^ [L 2 (u") + L 2 ( M ( 1 )) + 4L 2 ( 



+L 3 (u n y v n ) + L 3 (u ( - 1 \v^)+4L 3 (u^,v^)] 



(52) 



For |Aa 2 | S> 1, last inequality in (|47|) is the most relevant one. Firstly, we min- 
imize I — 1 + 3(1 — 2ci)(ci + 4c 2 )| (taking into account (|46[0 . Consequently, we 
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choose c 2 



4 ^3(1 - 2ci) 



ci I . This condition minimizes the factor accompa- 



nying A 3 a 2 in (|47[) . The remaining inequalities from (|46|) and (|47|) reduce to 



-17- v / 2377" -17+v / 2377 
- < Ci < - 



and 



Aa 2 -1 + 



72 



(-l + 6ci - 12c?) 



72 



(53) 



Secondly, we minimize 



3(2 Cl - 1) 

j-l + gci - 12c?) 



(taking into account its allowed range). 



3(2 Cl - 1) 
^ \/3 

The minimum is placed at ci = — - — . Therefore, the resulting values are 



(ci,c 2 ) 



3-V3 -l + \/3 
6 ' 8 



and the method is written as 



// ■- it™ + zi£ Li (u , v n ), 



v n + At 



(3 + V3) , „. (3-V3) r , ri)x , r / n m 
^ '-L 2 (u ) + ± '-L 2 {u K ') + L 3 (u ,v ) 



(55) 



(56) 



' = l[3u n +u^+AtL 1 (u^,v^)] 



(2) n i At 

v — v H — — 
4 



(^+^) L2(u «) + (zi+^) L2 ( u (i)) + | (3 - V / 3)L 2 ( W W) 



+L 3 (n",« n ) + L3^ (1) ^ (1) )] , 



n+l 



1 u n + 2u^ +2AtL 1 (u^,v^) 

A 

v n + ^ [i 2 (w n ) + L 2 ( W «) + 4L 2 (u 



(2h 



(57) 
(58) 



+L 3 (n", v n ) + L 3 (u^,v^) + 4L 3 (^ (2) ^ (2) )] 



As in the second-order PIRK methods, depending of the value of |Aa 2 | = 
|Aa 2 |Zli 2 , it will be more convenient to use a particular set of values for the coef- 
ficients. This fact will be illustrated in Sect.JU 



3.5 Comments on the PIRK methods and some IMEX ones 

In this subsection we comment briefly on the derived PIRK methods, the second- 
order IMEX-SSP2(2,2,2) and third-order IMEX-SSP3(4,3,3) from [5]. 

The IMEX-SSP2(2,2,2) method particularized to the system (TJ) reads: 

(1) _ n 

(1) ..n' , A+r t..n\ (59) 



V = V 



AtL 2 (u n ), 
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w (2) _ u n j_ M r> Ui n Jl) 



^ =v n + At [(1^2 7 )L 2 ( W ")+ 7 L 2 (^ 2 ))+L 3 K,« (1) )] , 1 j 



v n+1 = v n + 0.5 At 



L 2 (u n ) + L 2 {u^) + L 3 (u n ,v^) + L 3 (^ (2) V 2) )] , 

(61) 

where 7=1- 1/V2. The IMEX-SSP2(2,2,2) method is L-stable, and, in com- 
parison with the both derived second-order PIRK methods, there is an additional 
intermediate stage. On one hand, the second-order PIRK methods provide a stable 
evolution for a wave equation in spherical- type coordinates as it is shown in Sect.HJ 
On the other hand, the evolution of a wave-like equation with extra complex non- 
linear source terms using a numerical scheme with an extra stage is translated into 
extra computational costs. 

The IMEX-SSP3(4,3,3) method particularized to the system (JTJl reads: 



(i) 

(1) ? 



+ AtL 2 {u n ), 



(2) n 

(2) n 



(62) 

(63) 
(64) 



«W = u n + 0.25 At \l x (u n ,v n ) + Li(u w ,v w 

« (4) = v 11 + At \BL 2 {u n ) + (0.5 - a - (3)L 2 (u {3) ) + ci 2 (n (4) ) (65) 
+0.25 L 3 (u n ,v n ) + 0.25 L 3 (u (3) y 3) )] , 



« (3) =u n + AtL 1 (u n ,v n ), 

v {3) = v 11 + At \{1 - a)L 2 (u n ) + aL 2 (u {3) ) + L 3 (u n , v n ) 



n+l n | At 

6 



L 1 ( U ",7j™) + L 1 (^ (3) ,« (3) )+4L 1 ( lt (4) ,7jW)] 



n+l n . At [ . „. . (3), . , (4). (66) 

' = d + -g- \L 2 {u ) + L 2 (u^ ; ) + AL 2 (u y ; ) 



+L 3 («" , v n ) + L 3 (^ (3) , v {3) ) + 4L 3 (u^ , )] 



where a = 0.24169426078821, /3 = 0.18957643480295. Due to the particular struc- 
ture of the system ([TJ) , the resulting two first stages can be omitted, and the final 
number of stages is similar to the case of third-order PIRK methods. The trivial 
value for results from a cancellation, as a consequence of the presence of op- 
posite coefficients in the numerical scheme. Notice that negative coefficients are 
avoided in the derivation of the PIRK methods. 

The IMEX-SSP3(4,3,3) method corresponds to a third-order PIRK one, with 
C1 = a = 0.24169426078821 and c 2 = (1 - 3ci)/4 = 0.06872930440884. These 
values for the coefficients are close to one of the optimal sets deduced in Sect. 13.41 
ci = 4c 2 = 1/4. 
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4 Numerical experiments 

In this section, the PIRK methods described in Sect. [3] are applied to the case of 
the time evolution of a wave equation for a scalar, h, in spherical coordinates. A 
wave equation for h can be written as: 

d tt h = Ah, (67) 

where A is the Laplacian operator. Eq. (|67p can be rewritten as a first-order system 
in time, with the addition of an extra auxiliary variable, A, as follows: 

d t h = A, 

d t A = Ah. (68) 

In this case, according to system (fTJ), the variables can be identified as (u,v) = 
(h, A), and the operators as Ci(h, A) = A, C<2[h) = Ah and Cz{h, A) = 0. Therefore, 
a i = 7i = 72 = and c*2 = 1 in system The eigenvalues of the linearized 

explicit part are wi = uj2 = 1 and hence dex = 1 and trex = 2. A G R~ and its 
value depends on the particular solution for h and the discretization of the operator 

A. Spherical coordinates are used, being ( — -J^, — }— ^-J— ] the corresponding 



orthonormal basis. 

Eq. (|67[) has solutions of the form 



dr' r 88' r sin dtp 



h(r,e,tp,t) ~ ji{kr)Y lm (9,(p) cos kt, (69) 

being ji the spherical Bessel function of first kind of order I and Yj m the spheri- 
cal harmonics. The value of k, a positive real constant, is determined by imposing 
boundary conditions. We search for solutions inside a sphere of radius unity impos- 
ing h(r = 1, 8, ip, t) = 0. For fixed I and m it is possible to compute the eigenmode 
frequencies, k n [, as the zeros of the spherical Bessel function of order I, being n = 1 
the first zero and so on. 

We have performed ID-spherical, 2D-axisymmetric and 3D simulations of the 
system using as initial data solutions with n = 1 at t = 0. We use values of (l,m) 
equal to (0, 0), (2, 0) and (2, 2) for the ID, 2D and 3D case, respectively. In this way, 
the initial data is regular and fulfills the symmetries in each case. Furthermore, 
the data is non trivial for each symmetry in the sense that there are not trivial 
cancellations for 2D and 3D simulations. We consider symmetry with respect to 
the equatorial plane, 8 = ir/2. 

We use a finite difference scheme to solve the system using an equally-spaced 
grid with n r , ng and n v grid points in the r, 8 and tp directions, respectively. We use 
derivatives of cell-centered Lagrange interpolation polynomials [12] to compute the 
first and second spatial derivatives appearing in the Laplacian operator, achieving 
fourth, sixth and eigth discretization order. Boundary conditions are imposed by 
using a number of ghost cells consistent with the discretization stencil. At r = 1 
the analytical solution is imposed as boundary condition. At all other boundaries 
symmetry conditions are used. The maximum time step is determined by the CFL 
condition for the speed of the wave, which is 1. The time step in the simulations 
is a fraction smaller that the maximum time step, the CFL factor, in the interval 
[0,1]. 
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We can estimate the absolute error of the numerical evolution by comparing 
the numerical solution with the analytical one given by Eq. (|69|) . As a measure of 
the global error during the numerical evolution, we compute the Z/2-norm of the 
difference between the numerical and the analytical solutions at a given time, 



L a (A)(t)= Y [hnnm(r,e,(p,t) -hanaM,^)] 2 (fcr) 2 . (70) 

n r ngriu> \ * — ' 



4.1 Stability 

We have studied numerically the stability of the first, second and third-order PIRK 
methods. This study involves the numerical computation of a wide parameter 
space, including the coefficients c, of the PIRK methods and the CFL factor. Up 
to 10000 simulations have to be performed to cover this parameter space, so we 
have decided to use a single numerical setup as reference for the stability study. 
We use (n,l,m) = (1,2,0) for the initial data in 2D-axisymmetry with equatorial 
symmetry. In this case, k\2 ~ 5.763. We use (n r ,ng) = (100,32) grid points and a 
fourth-order spatial discretization scheme. 

With the numerical setup fixed, we can estimate the value of x := — A«2 = 
— \(At) 2 , which is the relevant quantity for the stability of the system. For solutions 
of the form (|69[l , Eq. (|68|) can be written as 

d t h = A, 

d t A = -k 2 h. (71) 

Therefore, A = -k 2 . The minimum value of k corresponds to the fundamental 
mode (n,l,m) = (1,0,0), i.e. fc m ; n = fcio = 7r. This sets the lower limit for x, 

x min = 7v 2 {At) 2 . (72) 

Note that in the limit of infinite resolution, the CFL restriction results in At — > 
and x min — > 0. We show in the next subsections that this is indeed the case for the 
typical resolutions used in practical applications. 

The maximum value of k corresponds to solutions with typical spatial varia- 
tions of the order of the smallest grid cell size, Al mln . For these solutions, 

femax ~ — rj \ i (^3) 

and hence the upper limit for x is 

z max = ^ ax (^) 2 w f ^^A (CFL factor) 2 = 4^ 2 (CFL factor) 2 , (74) 

being the CFL factor = At/Atmux and Atm&x = Al m i n the maximum time step al- 
lowed by the CFL condition, which coincides with the smallest grid cell size. Note 
that imu does not depend on the resolution. As a consequence, the stability prop- 
erties of the system do not depend on the resolution, but only on the CFL factor. 
For convenience we define x := x/(CFL factor) 2 ; therefore, x m i n = ir 2 (At maj x) 2 and 
imax = 4-7T 2 . By keeping the numerical setup fixed, we expect that the limits of x 
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stay constant in all our simulations. Effects of the dimensionality and resolution 
are discussed in section \5\ 

For each order of the PIRK methods, we compare stability behavior of the 
numerical simulations with the stability criterion of Sect. [3] for the determinant of 
the corresponding matrix of a given method. We compare this stability criterion 
with the general one for the eigenvalues of the system. 

4-1.1 First-order method 

The expression for det(Mi) particularized to the system (|68[) with x > 0, is 

det Mi = 1 + (1 - ci)x, (75) 
and the eigenvalues of Mi are 




The stability criterion, |e±| < 1, for x ^ 0, results in 

1 < Cl < | + |, (77) 
x < 4. (78) 

The condition | det Mi| < 1, used in Sect. [3j results in 

1 < ci < 1 + -, (79) 

x 

which is a necessary, but not sufficient, condition for stability. Conditions (|77[) - 
(|79[) are more restrictive for larger values of x, so the value x max determines the 
stability properties of the system. 

It is therefore relevant to study the stability properties of the numerical solution 
depending on the coefficient c\ and the time step At. We have performed series of 
simulations for several CFL factors (0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.85 and 0.9), varying 
the value of c\ from to 5 in steps of 0.001. We have evolved each simulation up to 
time t = 54, i.e. about 50 oscillations. The left panel of Fig. [1] shows the evolution 
of the Z/2-norm for a selection of simulations. The behavior of this quantity allows 
us to distinguish between numerically stable and unstable evolutions. Numerically 
stable evolutions (solid lines) show an oscillatory behavior of the Z/2-norm with 
values smaller than one, while numerically unstable evolutions (dashed lines) show 
an exponential grow of the I/2-norm which becomes much larger than 1. Hereafter, 
we define numerical stability for our numerical simulations as those in which L2- 
norm< 1 for a time corresponding to 50 oscillations, i.e. t = 1007r/w. 

For each CFL factor, we find a maximum and minimum values of c± such that 
all simulations within this two values are numerically stable, while lower or higher- 
values lead to numerically unstable evolutions. No numerically stable values of 
ci were found for a CFL = 0.9. The right panel of Fig. [T] shows the stability 
limits (red and black circles). We can compare these limits with the predictions 
of Eqs. (|77p and (|78p . For all CFL factors the lower limit is ci = 1 and coincides 
with the prediction of Eq. (|77|) . Following Eq. (|77[) . we fit the upper bound by a 
curve of the form pi + p2(At max / At) 2 , with p\ = 0.501, P2 = 0.3746 and At max 
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upper limit 
lower limit 




Fig. 1 Stability of the first-order PIRK method. Left panel: time evolution of the L2-norm 
for simulations with CFL factor 0.5 and ci values of 0.9 (blue), 0.99 (magenta), 1 (red), 1.5 
(orange), 2 (green) and 2.05 (black). Solid and dashed lines represent numerically stable and 
unstable simulations, respectively. Right panel: stability region depending on the values for 
ci and the CFL factor. Red and black circles are the maximum and minimum values of ci, 
respectively, for which a numerically stable evolution is found for each CFL factor. Solid lines 
are the boundaries of the stability region (orange area) given by Eqs. J77j and (1781 . for the 
fitted value of x = 5.340. The boundary of the region | det(Mi)| < 1, which partially coincides 
with the stability boundaries, is also plotted (dashed lines). 



the maximum time step predicted by the CFL condition. We confirm the I /(At) 2 
behavior of the upper limit and a value of pi compatible with the predicted i /2 in 
Eq. (|77p . From the fitted value of P2 we estimate x = 5.340, which is of the order of 
magnitude of i max ~ ty 2 . Using this estimation for x, it is possible to compute the 
maximum CFL factor using the condition (|78p . It results to be 0.8656, consistent 
with not finding any stable simulation for CFL = 0.9. 

If we consider the stability criterion used in section l3~2l (dashed-line in Fig.[l| 
for the determinant of Mi, we observe that this criterion is overestimating the 
stability region. However, the estimated optimal value, c\ = 1, lays inside the 
stability region and is indeed the value such that the maximum CFL factor is 
achievable. 

The first-order ERK can be recovered by setting c\ =0 0- Since none of our 
numerical simulations with c\ = and CFL factors between 0.3 and 0.95 show 
stable evolutions (see right panel of Fig. [1]), we have performed simulations de- 
creasing the CFL factor to try to find the stability limit. We were not able to find 
stable evolutions for CFL factors as low as 0.001. The reason for this behavior is 
that the eigenvalues of the first-order ERK, given by Eq. (|76|) with c\ = 0, read 

e± = 1 ± v 7 ^, (80) 

and the stability criterion, e±j < 1, is only fulfilled for x = 0, i.e. At = 0. If we 
decrease the CFL factor, i.e. approach je±j = 1, the instability appears at later 
time in the simulation (e.g., for CFL = 0.001 the instability appears at t ~ 7). The 

1 Alternatively, identical result can be obtained by setting ai = 71 = A = 0, C12 = 1 and 
71 ^ for all of the PIRK methods presented in this work. 
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consequence is that first-order ERK can be used for finite-time evolutions provided 
a sufficiently small CFL factor is used. However, the time step restriction of the 
first-order ERK is significantly larger (several orders of magnitude) than the one 
of the first-order PIRK. 



4-1-2 Second-order method 

The expression for det(M2) particularized to the system (|68[) with x > 0, is 

det(M 2 ) = i [4 - x{2c 2 - ci - 2cic 2 )] , (81) 
and the eigenvalues of M 2 are 
e± = l-f + |-(l-2e 2 )z 2 

± \\J-x [64 - 16(1 + 2ci - 2c 2 )x + 8ci(l - 2c 2 )x 2 - cf (1 - 2c 2 ) 2 x 3 }. (82) 
The stability condition, |e±| < 1, leads to 



- ( 1 - - ) <ci(l-2ca)< 

X \ X J X 

c± — c 2 < 2/x, 



c 2 (2ci-l)<- -1 + - 

X \ X 

< ci +2c 2 (ci - 1). 
The condition | det(M 2 )| < 1 is equivalent to 

< ci +2c 2 (ci - 1) < - 



(83) 
(84) 
(85) 
(86) 

(87) 



which coincides partially with the boundaries of the stability region. The condi- 
tions (|83p - (|86p and (|87p are more restrictive for larger values of x, therefore the 
value i mm determines the stability properties of the system. 

We have studied the numerical stability of the second-order PIRK by perform- 
ing simulations using ci S [—0.5,1.5] and c 2 £ [—0.5,1.5], varying the coefficients 
in steps of 0.02. We have used several CFL values (0.5, 0.7, 0.8 and 0.9). Fig. [5] 
shows the stability region on the (ci,c 2 ) plane. According to Eqs. (I83l) - (l86|) . the 
boundaries of the stability region depend only on the parameter x. In order to 
estimate x, we perform a x 2 -mhrimization of the difference between the numer- 
ically determined boundaries for all CFL factors and the theoretically predicted 
values. We obtain x = 5.373 with high significance (x 2 = 353.7 for 385 degrees 
of freedom). Note that, as expected, the value of x is very close to the one ob- 
tained in the first-order PIRK method. The boundaries for the fitted value of x 
(solid- lines in Fig. [2]) agree very closely with the numerically computed ones. As in 
the first-order PIRK, the criterion j det A/ 2 | < 1 overestimates the stability region 
(dashed- lines in Fig. [2]) . 

We can compare the numerical results for the stability region with the optimal 
values computed in Sect. 13.31 (ci,c 2 ) = (1/2,0) for \x\ <C 1 (black circle in Fig. [2]) 
and (ci,c 2 ) = \(2- V2,V2~ 1) for \x\ > 1 (star symbol in Fig. [2]). Both values 
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Fig. 2 Dependence of the numerically determined stability region (orange area) on the 
(ci,C2) coefficients for the evolution of h using a second-order PIRK method. The bound- 
aries (solid lines) of the stability region according to Eqs. (183 P — <T86l) for x = 5.373, and the 
boundaries (dashed lines) of the region |dct(Af2)| < 1, which partially coincide with the sta- 
bility boundaries, are plotted. The values for the optimal second-order PIRK method with 
( Cl ,c 2 ) = (1/2,0) (black circle) and (ci,c 2 ) = 1/2(2 - \/2, - 1) (star symbol), and the 
second-order ERK (black triangle) are also plotted. 



allow for stable numerical evolutions with CFL factors close to unity. According to 
the fit, x = 5.373 (CFL factor) 2 , which is larger than unity for CFL factors larger 
than 0.434. Indeed, the choice (ci, C2) = ^(2— \/2, \/ 7 2— 1) is better suited for higher 
CFL values, while the choice (ci,C2) = (1/2,0) becomes unstable. Therefore, we 
recommend the optimal values (ci, C2) = 2"(2 — a/2, V% — 1) as they seem to provide 
stable evolutions with the largest possible CFL factors. 

The second-order ERK correspond to the case (ci, eg) = (0, 1/2) and is uncon- 
ditionally unstable. In this case, the eigenvalues are 
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and the stability criterion, |e±| < 1, is only fulfilled for At = 0. Therefore, the 
observed numerical behavior of the second-order ERK is very similar to the first- 
order ERK. 

4-1.3 Third-order method 



The expression for det(Ms) particularized to the system (|68[) with x > is 



det(M 3 ) = 1 + — (ci - 4c 2 ) + _[_! + 3(1 - 2 Cl )( Cl + 4c 2 )], (89) 



and the eigenvalues of M3 are 

,.2 



e± = 1-1 + ^(1 + 01-402) 



24 



192 (i-3) 



-16x 2 (3ci(l - cl - 4c 2 ) + 1) + x 3 (l + ci - 4c 2 ) 2 ] 



1/2 



The stability criterion, [e±| < 1, for x ^ results in 



12 /-, 4\ , 12 

— 1- - < 1 + ci - 4c 2 < — , 

X \ X X 



1 + 3(2ci - l)(ci + 4c 2 ) < - ( — - 1 

a; \ a; 



l + 3(2ci-l)(ci+4c 2 ) < 



12 /4 



1 - 1 + 2(1 + ci - 4c 2 ) 



6 



(ci - 4c 2 ) < 1 + 3(2c a - l)(ci + 4c 2 ). 



(90) 

(91) 
(92) 
(93) 
(94) 



The condition | det(M3)j < 1 is equivalent to 

-(ci - 4c 2 ) < 1 + 3(2 Cl - l)(ci + 4c 2 ), 



-144 < 6x 2 (ci -4c 2 ) + z 3 [-l + 3(1- 2ci)(ci +4c 2 )], 



(95) 
(96) 



which coincides partially with the boundaries of the stability region. The condi- 
tions (|91|l - (|93p and (|96|) are more restrictive for larger values of x, therefore the 
value imax determines the stability properties of the system at these boundaries. 
However the condition (J94J) (which is equivalent to ([95])), can be more restrictive 
for small values of x depending on the value of c\ and c 2 . For this case, both i m ax 
and a; m in are relevant for the analysis of the stability. 

We have performed the same numerical stability analysis for the third-order 
P1RK as in Sect. 14.1.21 The numerically computed stability regions depending 
on (ci,c 2 ) are shown in Fig. [3] for different CFL factors. As in the second-order 
method, the stability regions shrinks for increasing CFL factors. We have per- 
formed a x 2 minimization to determine the value of x which fits better the bound- 
aries of the stability region with the expressions (|91[) -(|94 [ ). We obtain x = 5.322 
( x 2 = 109.0 for 255 de grees of freedom). For the numerical setup considered, 
we estimate x m i n = 6.945 x 10~ 7 <C 1, so we consider as a boundary the con- 
dition (|94[) setting x = (brown solid line in Fig. [3]). The resulting boundaries 
(solid lines in Fig. [3J fit very well with the numerically estimated stability region 
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c l c l 

Fig. 3 Dependence of the numerically determined stability region (orange area) on the (ci, C2) 
coefficients for the evolution of h using a third-order PIRK method. The boundaries (solid 
lines) of the stability region according to Eqs. ■ll' .l 1'iit for x = 5.322 (black lines) and x = 
(brown line), and the boundaries of the region |det(M3)| < 1, which partially coincide with 
the stability boundaries, are plotted. The values for the optimal third-order PIRK methods, 
(ci,c 2 ) = (1/4, 1/16) (black circle) and (ci,c 2 ) = ((3 - V3)/6, (-1 + V3)/8) (star symbol), 
and the third-order ERK (black triangle), are also plotted. 



(orange area in Fig. [3}. The condition |detMs| < 1 (dashed lines in Fig [3]) over- 
estimate the stability region; however, the optimal values obtained in Sect. 13.41 
(black circle and star symbol in Fig. [3]) lay inside the stability region for all CFL 
factors studied. Since both values are very close, any of them could be used in 
the case of the wave equation for CFL factors close to unity. If we compare to 
the third-order IMEX-SSP3(4,3,3) in [5 , corresponding to ci = 0.24169426078821 
and c 2 = (1 - 3ci)/4 = 0.06872930440884, the stability properties are very similar 
since the coefficients are very close to the optimal PIRK values, and therefore the 
method is numerically stable for the considered CFL factors. 
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Fig. 4 Numerical error estimation for series of ID (circles), 2D (triangles) and 3D (squares) 
simulations, normalized to the reference resolution of each series, as a function of the time 
step. Results for the first-order (black symbols), second-order (red symbols) and third-order 
(blue symbols) PIRK methods are shown. Black-solid line, red-dashed line and blue-dotted line 
correspond to power law fits for first-order, second-order and third-order PIRK respectively. 



The third-order ERK corresponds to (c±,C2) = (0,1/4) (black triangle in 
Fig- El- The eigenvalues of M3 in this particular case are 

e ± = l-|±iV-^-6) 2 - (97) 
The condition je±j < 1 results in 

< x < 3. (98) 

For the fitted value of x, this condition restricts the stability of the third-order 
ERK to CFL < 0.751. Accordingly, our numerical simulations show numerical 
stability for 0.5 and 0.7 CFL factors, but not for 0.8 and 0.9. We conclude that, 
although the third-order ERK is stable, the time step is more restricted as in the 
third-order PIRK. 



4.2 Convergence 

We have studied the convergence properties of the PIRK methods by performing 
series of lD-spherically symmetric, 2D-axisymmetric and 3D simulations. In each 
one, we use reference models with resolutions n r = 50, (n r ,ng) = (50,16) and 
(n r , rig, n v ) = (50,8,32), respectively. We increase the resolutions in all directions 
by the following factors: 2, 4, 8, 16 and 32 (ID); 2, 4, 8 (2D); 1.5, 2, 4 (3D). The 
I/2-norm of h after 1/4 period of oscillations, t = tt/u/2, is used as an estimation 
of the error. To ensure that we measure the error due to the time discretization 
and not due to the spatial one, we choose the spatial order carefully. The time 
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step in ID, 2D and 3D simulations in spherical coordinates scales as AtiD ~ Ar, 
At2D ~ ArAO, and At 3 D ~ ArAOAtp, respectively, i.e., if we increase the resolution 
in all directions by a factor 2, At decreases in a factor 2, 4 and 8, respectively. 
Therefore, to ensure that increasing resolution the spatial discretization errors 
decrease faster that the time discretization errors, we choose fourth-order spatial 
discretization for the first-order PIRK, fourth-order (ID) and sixth-order (2D, 3D) 
for the second-order PIRK, and fourth-order (ID) and eighth-order (2D, 3D) for 
the third-order PIRK. In all cases, CFL = 0.8. 

Fig. 3] shows the error, normalized to the error of the reference simulation of 
each series, as a function of the time step. Independently of the dimensionality 
of the simulation, the error falls with decreasing time step as expected from the 
convergence order of the PIRK used. Fitting to a power law, the convergence 
orders obtained are 0.87, 1.93 and 3.02 for first, second and third-order PIRK 
methods, respectively, very close to the expected ones. Note that for the 2D and 3D 
simulations using the third-order PIRK with highest resolution, the non-rescaled 
value of the Z/2-norm is close to 10 -14 and its value is limited by the machine 
accuracy. These two points are not converging anymore and we do not consider 
them to compute the convergence order. 



5 Conclusions 

In this work, PIRK methods, from first to third-order of convergence, have been 
derived to evolve in time systems of non-linear partial differential equations con- 
taining stiff source terms. The origin of the stiffness can be due to a particular 
choice of coordinates; in the case of spherical coordinates, 1/r and l/sin# factors 
appear in the equations and can be interpreted numerically as stiff terms close 
to r = and = 0,7r. Optimal SSP ERK methods are recovered when implicitly 
treated parts are neglected. The optimal values of the remaining coefficients are 
chosen according to stability criteria. Although the partially implicit parts confer 
stability to the system, no analytical or numerical inversion of any operator is re- 
quired, and the computational costs of the derived PIRK methods are comparable 
to those of ERK ones. 

We have applied the derived PIRK methods to the evolution of a scalar-wave 
equation and studied the numerical stability of the system depending on the co- 
efficients of the PIRK methods and the CFL factor. We were able to explain the 
boundaries of the stability region by using the bound on the eigenvalues of the 
system, |e±| < 1. The condition for the determinant, | det M\ < 1, which was used 
to optimize the coefficients of the methods in the first part of the present work, 
overestimates the size of the stability region. However, in all cases, the optimal 
values lay inside the stability region, and correspond to values for which the max- 
imum CFL factors can be achieved. Therefore, although results of the study of 
this particular system cannot be extrapolated to a general case, we think that this 
behavior could be similar for other systems of equations, and the optimal values 
obtained here can be use in more general applications with similar performance. 

The stability region of the methods can be determined in terms of the variable 
x. This value does not depend on the particular resolution of the system and can 
be used to check consistency of the derivation of the methods. We have determined 
the value of the upper limit of x for a set of 2D simulations with the same numerical 
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resolution and spatial discretization order. In all cases this value is very similar: 
5.340, 5.373 and 5.322 for first, second and third-order PIRK methods, respectively. 
These values are consistent with the growth of modes with spatial-scales of order 
the smallest grid cell size. In the third-order PIRK method, we were able to confirm 
that the lower limit of x is close to 0, as expected for the fundamental mode. 
Finally, we have checked the convergence order of our optimal PIRK methods, 
in agreement with the expected ones, in ID-spherical, 2D-axisymmetric and 3D 
simulations. 

Optimal SSP ERK methods can be viewed as particular cases for some values of 
the coefficients. First and second-order ERK result to be unconditionally unstable 
in both the linear stability analysis and the numerical simulations. The third- 
order ERK is stable, but the maximum CFL factor achievable is lower (0.751) 
than in the corresponding third-order PIRK (> 0.9). We suspect that the success 
of ERK methods in other works applied to wave- like equations could be due to the 
wide spread use of numerical-dissipation terms (e.g. |13jV which could avoid the 
growth of small-scale unstable modes. Our numerical method does not need these 
numerical-dissipation terms to get stable evolutions. We have observed similar 
performance in the comparison of the second and third-order PIRK methods with 
the IMEX-SSP2(2,2,2) and IMEX-SSP3(4,3,3) of [5], respectively. Regarding the 
second-order methods, the PIRK needs one stage less, which can be an advantage 
concerning the computational costs. Regarding the third-order methods, we have 
checked that the two first stages of the IMEX-SSP3(4,3,3) are redundant and this 
method corresponds to a PIRK one with coefficient values very close to the optimal 
ones derived in this work. 

These methods are also appropriate to evolve generalized complex wave equa- 
tions in spherical coordinates, as it has been shown in a recent application [14] for 
the evolution of Einstein equations in a particular formulation |15j using spherical 
orthonormal coordinates. The results of using these numerical schemes to evolve 
other systems of equations have to be explored, but this is beyond the scope of 
this work. 



A Stability conditions for second-order PIRK methods 

In this appendix we detail the derivation of the inequalities I I29II and (130 D . We would like to 
guarantee 

- 4 < Kx + Xa 2 K 2 {l - 2ci + 2c 2 ) + A 2 o|(2c 2 - ci - 2cic 2 ) < 4, (99) 

for u>i = 0,±1, i = 1,2, and in the cases IA02I 2> 1, Ac*2 ~ — ^ an< ^ l^ Q 2| "C 1: 

1. Upper bound, i^i = uj 2 = ±1: K\ = 4, K 2 = 0. In all the cases we require 2c 2 (l — ci) — ci < 
0. 

2. Upper bound, lui = —uj 2 = — ± 1: K\ = 4, K 2 = 2. 

For I Xa 2 | -C 1 , we require < 1 — 2ci + 2c 2 . This condition is sufficient for the cases 
|Act2| 2> 1 and \a 2 »s — 1. 

3. Upper bound, uj± = ui 2 = 0: K\ = 1, K 2 = 1. Previous derived conditions are sufficient for 
all the cases. 

4. Lower bound, uji = ui 2 = ±1: K\ = 4, K 2 = 0. In all the cases we require —8 < A 2 a|(2c2 — 
ci — 2ciC2). This condition will be satisfied by the following ones. 

5. Lower bound, u)\ = —u> 2 = — ± 1: K\ = 4, K 2 = 2. 

For |A«2| <SC 1, we require —4 < Ack2(1 — 2c+2c 2 ). 

For |A«2| S> 1, previous derived conditions are sufficient. 

For A«2 ~ _ 1, we require < 6 + 5ci — 6C2 + 2ciC2. 
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6. Lower bound, uii = u)2 = 0: K\ = 1, K2 = 1. 

For |A«2| "C L previous derived conditions are sufficient. 
For |Aq2| S> 1, we require —5 < A 2 a2(2c2 — c\ — 2c\C2). 
For Aa?2 ~ _ 1 1 we require < 4 + ci — 2ci c 2 • 



B Stability conditions for third-order PIRK methods 

In this appendix we detail the derivation of the inequalities II46H and 147 1 . We would like to 
guarantee 

_ x< Kl \a 2 K 2 {-l + c l -4c 2 ) 
~ 36 24 
A 2 a 2 

+ ^^[ci - 4c 2 + (dex - l)(4c 2 - c\ - 4cic 2 ))] 

-^[-1 + 3(1 -2ci)( C1 +4c 2 )]<l. (100) 

for uji = ±1, i = 1,2, and in the cases | Act2 1 S> 1, Ao?2 ~ —1 and |Aa 2 | <S 1: 

1. Upper bound, wi = UJ2 = 1: K\ = 36, K2 = 0, dex = 1. 

For |Aa2| <C 1, we require c\ — 4c2 < 0. 

For |Aa 2 | > 1, we require -1 + 3(1 - 2ci)(ci + 4c 2 ) < 0. 

For Aa2 ~ — 1, previous derived conditions in this appendix arc sufficient. 

2. Upper bound, lo\ = u>2 = —1: K± = 4, K2 = 8, dex = 1. 

For |A«2| li we require Aq?2(— 1 + c\ — 4C2) < 8/3. 
For |Aq2| S> 1, previous derived conditions are sufficient. 

For Ac*2 ~ —1, we require < 41 + 18(ci — 4c2) — 3(1 — 2ci)(ci + 4c2). Taking into 
account previous derived conditions, it is sufficient to require —20/9 < c\ — 4c2- 

3. Upper bound, w\ = 1 = — u>2 or oj2 = 1 = —uif. K\ = —12, K2 = 8, dex = —1. 

For |Act2| <SC 1, previous derived conditions are sufficient. 
For |A«2| 1, previous derived conditions are sufficient. 
For Aa 2 ft: -1, we require < 73 + 18c 2 - 180c 2 + 9ci(3 + 8c 2 ). 

4. Lower bound, 1^1=^2 = 1: K\ = 36, K2 = 0, dex = 1. 

For |Aq2| < 1, we require — 24 < A 2 a|(ci — 4c2). 

For |Aa 2 | > 1, we require A 3 a|[-1 + 3(1 - 2ci)(ci +4c 2 )] < 144. This condition will 
be satisfied by following ones. 

For Aa 2 ~ -1, we require < 9ci - 12c 2 - 6c 2 - 24cic 2 + 143. 

5. Lower bound, u)i = UJ2 = —1: K\ = 4, K2 = 8, dex = 1. 

For |Aq2| <SC 1, previous derived conditions are sufficient. 

For |Aa 2 | > 1, wc require A 3 a|[-1 + 3(1 - 2ci)(ci + 4c 2 )] < 80. This condition will 
be satisfied by following ones. 

For Aa 2 ~ -1, we require < 103 - 15ci - 6c 2 + 84c 2 - 24cic 2 . 

6. Lower bound, uii = 1 = — u)2 or ui2 = 1 = —uii: K\ = —12, K2 = 8, dex = —1. 

For |A«2| <S 1, previous derived conditions are sufficient. 

For |Aa 2 | > 1, we require A 3 a|[-1 + 3(1 - 2ci)(c! +4c 2 )] < 48. 

For Ac*2 ~ —1, we require < 6c 2 — 15ci + 36c2 + 24ciC2 + 71. 
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